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The detection and estimation of gravitational wave burst signals, with a priori unknown polar- 
ization waveforms, requires the use of data from a network of detectors. For determining how the 
data from such a network should be combined, approaches based on the maximum likelihood prin- 
ciple have proven to be useful. The most straightforward among these uses the global maximum 
of the likelihood over the space of all waveforms as both the detection statistic and signal estima- 
tor. However, in the case of burst signals, a physically counterintuitive situation results: for two 
aligned detectors the statistic includes the cross- correlation of the detector outputs, as expected, 
but this term disappears even for an infinitesimal misalignment. This two detector paradox arises 
from the inclusion of improbable waveforms in the solution space of maximization. Such waveforms 
produce widely different responses in detectors that are closely aligned. We show that by penalizing 
waveforms that exhibit large signal-to-noise ratio (snr) variability, as the corresponding source is 
moved on the sky, a physically motivated restriction is obtained that (i) resolves the two detector 
paradox and (ii) leads to a better performing statistic than the global maximum of the likelihood. 
Waveforms with high snr variability turn out to be precisely the ones that are improbable in the 
sense mentioned above. The coherent network analysis method thus obtained can be applied to any 
network, irrespective of the number or the mutual alignment of detectors. 



I. INTRODUCTION 

Several interferometric gravitational wave (GW) de- 
tectors are now operational around the world 0, H| 
in addition to an already operating set of resonant mass 
detectors @ . Since gravitational wave detectors can mea- 
sure both the amplitude and phase of incoming signals, it 
is possible to combine the data from such a network of de- 
tectors to not only detect GW sources better, compared 
to single detectors, but also estimate the waveforms of 
the two independent polarization components of the GW 
signal along with the direction to the source. 

A straightforward approach to signal detection and 
estimation using data from a network of detectors is 
through the likelihood functional. The likelihood 
functional of data x is the joint probability density func- 
tion p(x|0), where is a vector of parameters whose true 
value for any given instance of x is not known a pri- 
ori. In the case of GW detection, x is the collection 
of outputs from all the detectors in a network and the 
parameters are the sky position of the GW source and 
the waveforms of the two GW polarization components, 
h+{t) and hx{t). The problem of signal detection is con- 
cerned with deciding whether x is a realization of p(x|0o) 
or p{x\6 ^ 9o), where 6*0 denotes the signal absent case. 
The estimation problem deals with infering the true value 
of 6 once a detection has been made. According to the 
maximum likelihood approach to inference, the value 
^max at which the likelihood functional attains its max- 
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imum serves as an estimate of the true value of 6. The 
maximum, over 9, of the likelihood ratio (LR) functional 
p{x\9)/p{k\9o) yields a detection statistic. Thus, the like- 
lihood based approach conveniently yields both a detec- 
tion statistic and source parameter estimates. 

For binary inspiral signals, the polarization waveforms 
can be computed in terms of a small set of parame- 
ters describing the binary, such as the masses of the 
two stars, the orbital inclination etc. Likelihood based 
network analysis is fairly well understood in the case of 
such known waveforms 0, Q- In this paper, we con- 
sider likelihood based network analysis for burst signals 
that are short duration deterministic signals but have 
unreliable or simply unknown waveforms. Astrophysical 
sources of burst signals include Supernovae |^, Gamma 
Ray Bursts 10] and the merger of compact objects such 
as Neutron stars and Black Holes [Tl| . 

Although we do not have prior knowledge of burst sig- 
nal waveforms, it is possible to extend likelihood based 
network analysis to burst signals in a formal sense. The 
basic idea [iJiliSllJillS is to treat each sample of a burst 
signal as a parameter and maximize the likelihood over 
all the parameters. We refer to the resulting detection 
method as standard likelihood (SL). Other approaches 
to network analysis for burst signals have also been pro- 
posed in the literature 'T^,'l?|. The method proposed 
by Giirsel and Tinto is based on extremizing a functional 
that is constructed purely out of the response of each de- 
tector in a network and does not explicitly depend on the 
polarization waveforms. A minimum of three misaligned 
detectors are required in order to use the Giirsel-Tinto 
method whereas likelihood based methods do not have 
this limitation. 

Although it is quite elegant, the SL method has a fun- 
damental problem that was first noted in and system- 
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atically analyzed in [TBI (henceforth, KMRM). Called the 
two detector paradox, the problem is most easily stated 
for a network of two detectors. When the detectors are 
perfectly aligned, the SL detection statistic naturally in- 
cludes the cross-correlation of detector outputs. Since 
the response to a given GW source is identical in aligned 
detectors, this term serves to increase the sensitivity of 
the overall detection statistic. However, it disappears 
when the SL approach is used for even infinitesimally 
misaligned detectors. This is counterintuitive since we 
expect the cross-correlation term to remain important 
for small misalignments. 

The origin of the two detector paradox is that all wave- 
forms for h^{t) and h^(t) are allowed as solutions when 
maximizing the likelihood even though some of them 
would be improbable in that they will produce widely 
different responses in two nearly aligned detectors. It fol- 
lows that this problem can be solved by removing such 
improbable waveforms from the space of allowed solu- 
tions. In KMRM, this was achieved by imposing certain 
constraints on the space of waveforms yielding the class 
of constraint likelihood methods. Modifying the likeli- 
hood functional with a waveform dependent weight fac- 
tor, such as a Bayesian prior |l8j| . is another way of sup- 
pressing undesirable solutions. In the following we use 
the term "restriction" to encompass all approaches that 
remove improbable solutions or reduce their influence. 

The non-trivial aspect of setting up restrictions is that 
the performance of the resulting statistic must not be de- 
graded and, if possible, should be better than the SL ap- 
proach. Additionally, the restrictions must be applicable 
to all astrophysical burst sources otherwise the resulting 
statistic will be better for a particular type of source but 
may be worse for others. If sufficiently reliable knowl- 
edge for a p articular class is available, then it can be 
used |20| as an additional set of restrictions to further 
enhance sensitivity. The constraints proposed in KMRM 
(see Section Un)) satisfy the above condition, namely inde- 
pendence from assumptions about source dynamics. Sim- 
ulations show that constraint likelihood methods perform 
significantly better than SL, especially in reconstructing 
the sky position of a burst source. However, the con- 
straints used were a first example and are not necessarily 
unique or optimal. It is important to investigate indepen- 
dent avenues of formulating generic restrictions in order 
to gain more insight into this promising approach to net- 
work analysis. 

We investigate the issue of generic restrictions from a 
physical point of view. It is found that the more im- 
probable a set of waveforms /i+ and hx, the faster is 
the change in detectability of the corresponding source 
as a function of sky position. By measuring the rate of 
change of signal-to-noise ratio (SNR) as a function of sky 
position, therefore, we can restrict solutions that are un- 
desireable in the maximization of the LR. Monte Carlo 
simulations show that the detection statistic that results 
from the above restriction performs roughly the same as 
constraint likelihood methods. However, there exist im- 



portant qualitative differences in how the two methods 
process data. For example, unlike constraint likelihood, 
the divergence from SL of the method introduced here 
depends on the actual data. 

The paper is organized as follows. We review the stan- 
dard likelihood method in Section ^1 establishing our 
basic notation and conventions on the way. Section IIIII 
contains a discussion of the two detector paradox. Sec- 
tion ll V I considers the SNR variability of signals, discusses 
how it affects their detectability and defines a measure 
for it. Section contains the derivation of the detec- 
tion statistic and source parameter estimators that re- 
sult from incorporating SNR variability. Results from 
simulations are presented in Section IVll 



II. STANDARD LIKELIHOOD ANALYSIS FOR 
BURST SIGNALS 

We fix a Cartesian coordinate system with its origin 
at the center of the Earth, Z axis pointing at the North 
Pole and some arbitrary convention for the X and Y 
axes. This will be our fiducial frame of reference. Let the 
sky position of an incoming GW signal in this frame be 
denoted by 9, the polar angle, and 0, the azimuthal angle. 
An incoming GW signal is most conveniently expressed in 
the transverse traceless (TT) guage 21] associated with 
the source direction. The Z axis of the TT frame is 
oriented along the direction to the GW source. We adopt 
the convention that the X axis of the TT gauge lies in the 
plane formed by the TT frame Z axis and the Z axis of 
the fiducial frame. The GW signal is specified by its two 
independent polarization components /i+(i) and hx{t) in 
the TT frame. 

In this paper, we take into account the different orien- 
tations and geographic locations of currently operating 
interferometric detectors but treat them as identical in 
sensitivity. For instance, a GEO detector in the paper 
will mean a detector having the same location and orien- 
tation as GEO600 but with sensitivity similar to that of 
LIGO detectors. This simplifies the algebra considerably 
whithout altering the main result which is the relative 
comparison of methods. 

Since real data is digital, we shift to the discrete 
time domain and from here on denote any time depen- 
dent quantity g{t) as g[j], where t = jS, 5 being the 
sampling interval. A finite sequence of time samples 
{^[0], g[l], . . . , g[m — 1]} will be denoted by g. The scalar 
response of a detector s[k] to an incoming GW is given 
by, 

s[k]^W'^[k]D,j , (1) 

where W^^ is the symmetric trace free tensor that de- 
scribes the GW in the transverse-traceless gauge and -D^- 
is an equivalent tensor associated with the detector |2a |. 
(The Einstein summation convention is in force in the 
above expression.) The response Si[k] of the i"^ detector. 
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1, . . . ,Nd, is expressed as, 



s^[k]=K 



)h+[kd-S,]+F^,,i0,^)h^[kS-S,] , (2) 



where F^^i and Fx.i depend on the rotations involved in 
going from the TT fr-ame associated with a detector to 
the wave fr'ame associated with the fiducial frame and Si 
is the travel time for the signal from the origin of the 
fiducial frame to the detector. For a given 9, </>, we can 
treat the detectors as co-located {Si = 0) at the origin 
of the fiducial frame since the output from each detector 
can be shifted in time to compensate for Si. 

For a network of detectors, it is convenient to define 
quantities that are analogous to those for a single de- 
tector. Thus, we can define a network response vector, 



S[fc] = 



si[k] 



sjVd [k] 



(3) 



Vectors in M.^'', the space of real A^^-tuples, will be de- 
noted by boldface capital letters. By convention, ev- 
ery vector A G U^'^ will be a column vector. We de- 
fine the Euclidean scalar product in this space, A-B — 
A^B. A sequence of vectors will be denoted by A = 
{A[0],...,A[m-l]}. 



For given 6, 



S[k] = T{0,i 



h+[k] 

hy. [k] 



(4) 



where the matrix J- is defined as, 

/ F+.(9,(t>) Fx,i( 



(5) 



Given a network response vector S[fc], the polarization 
components can be obtained as. 



h+[k] 



(6) 
(7) 



provided 7i~^ exists. (Note that H is a symmetric 2-by-2 
matrix and, hence, has two real eignevalues.) 

Let the noise in the i^^ detector be n^. In the following 
we will assume that ni is a zero mean, white, Gaussian 
process, i.e., E[n4j]] = and E[n4j]n„J/c]] = SimSjk- In 
analogy to the network response vector S[A:], we define a 
network data vector X[A:], 



(8) 





The natural logarithm of the likelihood ratio (LR), 
A(X|6', (p, S), for an Nd detector network is given by 

, NaN-l 

A(X|0,0,S) = -^^(x.bf -(a;.b-]-s4j]f)(9) 

i=l j=0 
N-l 

= E 

J=0 



X[j].Sb-]-i||Sb-]f 



(10) 



where the detectors outputs are implicitly time shifted 
as explained above. 

The maximum likelihood ratio approach can be ex- 
tended to burst signals by substituting Eq. 01 into Eg .1101 
and maximizing A(X|0, 0, S) over each hx[j] in- 

dependently. Since the resulting equations involve only 
quantities at the same sample j, this involves maximizing 
the summand in Eq. 1101 independently for each j . From 
now onwards we will focus on only one time instant and 
drop the time index when there is no scope for confusion. 
The quantity to be maximized is, 



Ab-]=Xb-].S[j]-i||Sb-] 



(11) 



over the vector S in R^"*, the space of real A'^^-tuples. Let 
the maximum of A occur at Smax ■ If the maximization 
allows all possible S, then it is obvious that Smax = X. 
However, we are not allowed to choose any arbitrary S 
in the evaluation of A since from Eq. ^ S is a vector 
that depends on only two parameters and hx[j]. 

This implies that S must lie in a two dimensional plane 
determined completely by the antenna pattern functions 
of the detectors and passing through the origin of M.^"^ . 

Though the maximization of A is restricted to S lying 
in the response plane, in the maximum likelihood ratio 
approach there are no further constraints on the choice 
of S within this plane. We call the method of maximiz- 
ing A over S with no further restrictions as the standard 
likelihood (SL) method [H [H El El El . We denote 
the value of A (Eq. 110(1 , obtained by summing the max- 
imized values of A corresponding to each time instant, 
by Asl(0, 0). This is the detection statistic of the SL 
method. 



III. THE TWO DETECTOR PARADOX 

In the two detector case, the network data vector X lies 
in the response plane. When the detectors are aligned, 
we are constrained to use a network response vector S 
lying along the 45° line since the responses of the two 
detectors must be identical, si[k] = S2[k] — s[k] for all k. 
Maximizing A over the response s[k] yields. 



1/ 



A CX -(Xi +X2)'^ , 



(12) 
(13) 
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where A denotes the maximum value and A = A(s'). It 
follows that, 

N-l 

Asl{0,4>) a (^ibf +2^2bf + 2a;ib1x2b1) . (14) 

j=o 

Thus, the SL detection statistic contains a cross- 
correlation term, Xi[j]x2[j], which agrees with our 
intuition: since the responses of the two detectors are 
identical, the cross-correlation term will always have a 
positive definite mean in the presence of a signal and, 
therefore, should be useful for detection. 

However, the same procedure, when carried out for 
even infinitesimally misaligned detectors, yields Ag^ that 
has no cross-correlation term, 

si = xi ; S2 = X2 , (15) 
A cx x\+ x\ . (16) 

This clearly violates what one expects on physical 
grounds since the cross-correlation term should remain a 
useful discriminant against noise when the detectors are 
only infinitesimally misaligned. However, for the class of 
probable waveforms, the SL statistic becomes less sensi- 
tive than one in which a cross-correlation term is added 
by hand. Thus, the SL approach fails to deliver a phys- 
ically acceptable detection statistic in the two detector 
case. Following KMRM we call this problem the two 
detector paradox. 

The origin of the two detector paradox lies in the max- 
imization of A over all waveforms, and h^, including 
those that can produce very different responses si and S2 
for even infinitesimally misaligned detectors. Such wave- 
forms can lead to small or even zero cross-correlation 
between the responses. Physical intuition suggests that 
such waveforms should be imrobable. Mathematically, 
however, such waveforms are allowed since the matrix 
T is invertible, even if it is nearly singular, and some 
h+ and hy, can always be obtained for any given si and 
S2. As a result, the SL approach "throws" out the cross- 
correlation term since this term can now contribute pure 
noise. Geometrically, the two detector paradox comes 
about because the data vector X lies in the response 
plane and the SL approach does not constrain the choice 
of S in this plane. Hence, the solution S = X that re- 
moves cross-correlation is an allowed one in this case. 
This suggests that restricting the allowed solutions in 
the response plane can resolve the two detector para- 
dox. The SL detection statistic recovers pair-wise cross- 
correlations in the case of three or more misaligned de- 
tectors since X now lies outside the response plane. How- 
ever, the two detector paradox is just an extreme mani- 
festation of a problem that occurs for any arbitrary net- 
work: depending on the geometry of the network, some 
solutions are improbable and their inclusion in the max- 
imization of the LR reduces its sensitivity. 



IV. VARIABILITY OF SNR OVER THE SKY 

Though our physical intuition suggests that an in- 
finitesimal misalignment between two detectors should 
not result in very different responses si ~ S2, this does 
not forbid the existence of such sources. As discussed 
above, mathematically, h^{t) and hx{t) can be com- 
pletely arbitrary and, in particular, can be such that si is 
very different from S2- However, we expect such sources 
to be improbable. Is their a quantitative basis for our 
expectation? In this section we show that, indeed, there 
is a natural argument that "probable" GW sources, for 
which si ~ S2, should occur more often than "improb- 
able" ones. Note that the categorization of a source as 
probable or improbable also depends on the mutual align- 
ment of detectors. For example, two detectors oriented at 
45° to each other can often have very different responses 
for the same source. As shown below, our argument takes 
this factor into account. We first present the two detector 
case and then consider an arbitrary network. 



A. Two detector case 

Consider a gravitational wave source located at {9, cj)) 
with signal polarization components /i+[A:] and hx[k]. 
Let the detector response vector be S[k]. Now consider 
an identical source located at {9',4)'). By this we mean 
that the fiducial frame, along with the GW source, be ro- 
tated rigidly leaving the detectors in their original place. 
This leads to a network response S'[k] that is related to 
S[k] by, 

s' ^ T{e',(i)')T-\e,(t))s . (17) 

Now, assume that we use matched filtering l^sj to detect 
the source at the instant k. Thus, our detection statistic 
will be X[fc].S[fc] and X[fc].S'[fc] for the two sources re- 
spectively. The detectability of a source under matched 
filtering is measured by the signal to noise ratio (SNR), 
p, defined as the ratio of the mean of the detection statis- 
tic to its standard deviation. For the two sources above, 
the respective SNR wih be p = ||S[A:]|| and p' = ||S'[fc]||. 
Note that we are considering a more restricted version of 
matched filtering than the conventional one in which the 
detection statistic would be ^X[A:].S[A:] and the SNR 
would be = ^ ||S[A;]|p. We return to this point later 
in this section. 

We denote the relative difference in SNR for the two 
source locations as, 

e(0,0,0',0';^) = ^— ^, (18) 
P 

where S is the unit vector oriented along S. For small 
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displacements, 6' = 9 + 59 and 4>' 



we get, 



de 



(19) 



de_ 



d_ 

89' 



-,jm9\<j,')T-^i9,^)S\\^^ 



S 



f dT{9 



\ del) 



'-^T-^{9,<p)S] . (20) 



The norm of the gradient 



1/2 



(21) 



Te — dT/d9 and ~ dJ- /dcj), is a convenient measure 
of the variabiUty of the SNR of a source as a function of 
its location. ^ 

A high value of Q{9, 0; S) implies that the source would 
have to occur at very special positions on the sky to be 
detectable. This is because if it has an SNR above some 
detection threshold at one location, it may rapidly fall 
below it if it were to occur at a slightly different loca- 
tion. Consider two classes of sources, class A with a 
large associated Q{9, 0; S) and class B with low values 
of the same. Assume that events from both classes oc- 
cur with the same comoving rate and that they have the 
same distribution on the sky. Then, due to the rapid 
fluctuation of the SNR for source class A, the area of the 
sky over which its members will be detectable would be 
significantly less than that of class B. It follows that the 
observed rate of the class A sources would be much re- 
duced in comparison to class B. Therefore, even though 
the two sources may be equivalent in terms of, say, the 
energy emitted in GWs, the high SNR variability source 
would be detectable less often. Since a larger volume 
of space would required to detect such sources, it will 
naturally be harder to detect at a given false alarm rate. 

We thus arrive at the main point of this section. Given 
a pair of misaligned detectors, we should include only 
those responses S in our search for the maximum of the 
LR for which Q{9, </>; S) is low in some well defined sense. 
In this way, the detectability of sources with an intrin- 
sically low observable rate will be reduced but that of 
more frequently observable sources will be enhanced. It 
is important to note that the above argument does not 
require any prior assumptions about the astrophysics of 
sources. Also note that Q depends only on the orientation 
of S and not its magnitude. That is, it is only concerned 
with the shape of the waveforms of /i+ and /ix and not 
the distance to the source. 

Figure n shows polar plots of G{9,(j);S) as a function 
of S for fixed 9, (j) and different detector alignments. It 



is assumed that the detectors are colocated at the origin 
of the fiducial frame and coplanar with angles between 
their respective X arms of a. The antenna patterns are, 
in this case. 



^(1 + cos^ ^?) cos2(/) COS0 sin2(/)^ 

i(l -I- cos^ 6*) cos 2(0 -I- a) cos6i sin 2(0 -I- a) 



(22) 

Each plot is over all orientations of the detector response 
vector S. First, we see that for almost coaligned detec- 
tors, a = 0.1°, If we limit the search for the maximum of 
A (c.f., Eq.ini) over S such that 0(9, 0; S) < 10, say, then 
only detector responses with orientation very close to the 
45° and —45° line will be allowed. This means that the 
instantaneous response in the two detectors should be 
nearly equal. From this example, we see that waveforms 
with low SNR variability are also the ones that lead to 
nearly identical responses in the two detectors and, as dis- 
cussed above, are the ones which have relatively higher 
chances of being detected. This bears out our physi- 
cal intuition that sources that produce very different re- 
sponses in nearly aligned detectors should be rare. As a 
is increased, one sees that not only does G decrease but 
that it also approaches the same value for all orientation 
and, hence, all orientations of S are now allowed. The 
pattern opens up fully when the detectors are at 45° to 
each other. As a approaches 90°, the responses are again 
increasingly restricted around the ±45° lines. 

Though Fig. n shows that SNR variability is a natural 
measure for restricting the orientation of 5*, the particular 
measure G{9, (j>; S) allows both ±45° lines as the locii of 
restriction. For nearly aligned detectors, however, we ex- 
pect detector responses to be in-phase which means that 
the responses should lie in the vicinity of the +45° line. 
We introduce another function, T{9,(j); S), that is also a 
measure of SNR variability but has better properties in 
this respect. 



r2(0,0;5) 



(23) 



In fact, from the Cauchy Schwartz inequality, it follows 
that 



g{9,^;S)<Ti9,4>;S) 



(24) 



Hence, T{9, 0; S) is an upper bound on SNR variabil- 
ity. Figure 121 shows polar plots of T{9,(j); S) for fixed 9, 
(j) and different misalignment angles a. Now, for small 
misalignement (a = 0.1°), one sees that only responses 
along the +45° line are allowed while for misalignement 
close to 90°, only responses that are close to —45° are 
allowed. It makes sense, therefore, to use T{9, (p; S) as a 
measure of SNR variability instead of Q{9, 0; S). 

In the above discussion, we considered the detection of 
a signal at a single instant using a matched filter, leading 
to the detection statistic X[fc] ■ S [k] . Burst signals can last 
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FIG. 1: Plots of Q fEg. I2H for different misalignment angles. 



over several sampling intervals and, in general, the detec- 
tion statistic would in fact be a sum -^[^1 ' over a 
range oi k. In this case, a high SNR variability response 
S[/c] can be balanced by a very low variability response 
S[/c'] at some other instant, thus leading to a low variabil- 
ity for the overall SNR = E ||Sf ]^/^ This would allow 
a broader range of signals to be used for maximizing the 
overall likehhood A (c.f.. Eg. [Tn|) . Our choice of using 
instantaneous SNR variability is simply a conservative 
one since signals with low instantaneous SNR variability 
will necessarily have low variability of their overall SNR. 
Moreover, as discussed earlier, physical intuition suggests 
that the responses in closely aligned detectors should be 
nearly identical at each instant of time and not just on 
the average. 



B. General network of detectors 

We will now generalize the SNR variability measure 
r(6', S"), introduced in Eq. ESI to the case of an arbi- 
trary network of detectors. The relative difference in the 
observed SNR of the same source between two positions 



{9,(1)) and [e',cj)') is, 

e{B, 0, B', S) - \\T{e', ^')H~\e, <j))J'^{e, <j))S\\ - 1 . 

(25) 

As before, for a small displacement of the source. 



09' 
de 



(26) 
(27) 



and our measure of SNR variability will be defined as. 



r2(^,</,; S) 



TeU-^T^ S 



J'd>'H T S 



(28) 



Here, S" is a unit vector in the R^'' space but lying in the 
response plane. 

For algebraic simplicity, it is convenient to express the 
response vector S in terms of a basis in the response plane 
rather than the full i?^"* space. In order to do so, we need 
to fix two orthonormal vectors % and 5i on this plane. 
If Co and ei are the two orthonormal eigenvectors of 7i, 
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FIG. 2: Plots of r ffig. I23II for different misalignment angles. 



corresponding to eigenvalues and vi respectively, then 
it is easily seen that for, 

90 = -^HO,^)eo (29) 

V ko 

91 = -^Ti0,cj))ei (30) 
qo-qi = e^Hei = . (31) 

ll?o|| = l|gi|| = i. (32) 

Since, qq and qi have the same form as the detector re- 
sponses (c.f. Eq. they are confined to the response 
plane. Let 7 be the angle between S and qo. Then, 

S = ||S|| ($ocos7 + qisin7) , (33) 

The SNR variability measure T{d, ip; S) then becomes a 
function r{9, 0; 7) of 7. 

Recall Eq. ^2 for the definition of the quantity A that 
we want to maximize over the response vector S. We can 
now rewrite it as, 

A = (X • go) /OCOS7 + (X • $1) psin7 - , (34) 



where p = ||S||. From the above it is obvious that the 
solution for S that maximizes A is the projection of X 
onto the response plane. 



V. PENALIZED MAXIMIZATION OF 
LIKELIHOOD RATIO 



We have shown that SNR variability, as measured by 
r{6, (j); 7) defined in Eq.[2Hl correlates well with how im- 
probable a solution of the LR maximization problem is. 
Improbable waveforms /i+[fc] and [k] that lead to very 
different responses in nearly aligned detectors also corre- 
spond to a fast variation in the detectability of the source 
as it is displaced on the sky. Such sources naturally be- 
come rare compared to probable ones. In this Section, we 
construct a detection statistic that incorporates T(6, (j)] 7) 
to restrict the influence of improbable solutions on the 
maximization of the LR. 
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A. Algebraic simplifications 



fhi respectively. From Eg. 1341 and 1421 it follows that, 



We begin by making some algebraic simplifications. 
Substituting Eq. 123 into |2H1 we get, 



r2(0,0,7) = (cos7 sm7)A^(3°^';^) 



M 

A 
B 
C 



A C 
C B 



^(||.Feeo|p+||J-^eo|P) 



\\^4>ei\ 



(35) 
(36) 
(37) 
(38) 

(39) 



We now perform another change of basis in the response 
plane. Let the eigenvalues and eigenvectors of A4 be /ig. 
Ml ^ Mo ^nd mo and fhi respectively. Since fho and rhi 
form an orthonormal basis, 



cos 7 
sin 7 



fho cos ip + fhi sin ij; 



(40) 



for some angle tp between the response vector S and fho. 
Then, 



T^{6, (/); 7) = Mo cos^ + Mi sin^ tp 



(41) 



Henceforth, we denote r^{6, 0; 7) as T'^{0, 4>] ijj). 

Using the orthonormal basis vectors fho and mi, we 
can state 



X • $0 
X • qi 



a^fho + Pxfhi 



(42) 



where ax and Px denote the components along fho and 



p(cos7 sin7) (^x^|°) -^m' 
p (ax cosip + Px sin-0) - , 

plx cos {lp -(f)x) -^P^ , 



Ix = y/al + (31 , 



tan- 

V ax 



(43) 
(44) 
(45) 



The maximization of A has to be performed over p and 
-0. Thus, the SL solution \s p — Ix, the length of the 
projection of X onto the response plane, and ip = 4'x 
which is the angle between the projection and the fho 
vector. 

B. Penalized maximization 

Now we introduce one possible way to incorporate 
T{9,<pi,ip) into the maximization of A. We call this 
method penalized maximum likelihood ratio (PMLR). We 
propose to maximize A', 



)^ = \~aV^{e,(P,iP) 



(46) 



where a > is a parameter under our control. The 
extra term aV'^{9,(p,'ijj) acts as a penalty that prevents 
the solution from readily converging to the data vector X. 
This resolves the two detector paradox (c.f.. Section llTT|l . 

Once the solution p, ip that maximizes A' in Eq. 1461 is 
found, we find the value of A at this solution (not A'). 
Our final data functional will be the sum of the values 
of A, found as above, over all samples of the data X. 
Reintroducing the time index and denoting the resulting 
functional as Apl(0, 0), we get 



J 



A 



E 



max 



plx M] C0s(V'[j] 



y{Po cos^ tp[j] + Ml sin^ V'[i]) 



(47) 



r 



The global maximum of Apl(^, </>), over 9 and (p, will 
serve as the final detection statistic. 

The solution p, ip that maximizes A' is easily found. 



ta,n{2ip) 



p = lxCos{tp - (px) , 
sa\{2(px 



cos{2(px) - 2a{po - Mi)/^i 



(48) 
(49) 



Note that when a = 0, we recover the standard likelihood 
ratio statistic Ksl{S, (p) wherein the solution for S is the 



data X projected on the 2D response plane. Also, the 
importance of the penalty term diminishes as /^T^. If a 
signal is present, the latter increases with an increase 
in the instantaneous signal amplitude. Hence, whenever 
the signal amplitude is strong enough, the effect of the 
penalty term is reduced at that instant and the wave form 
estimators ip and p converge to the SL estimates. 

The term po ~ Mi: which is always positive by defini- 
tion, is a measure of how misaligned the detectors in a 
network are. This difference is reduced with an increase 
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in misalignment, hence reducing the effect of the penalty 
term. Fig. O shows hq — Hi as a function of sky position 
for two networks: (1) three LIGO detectors: one at Liv- 
ingston, and two at Hanford, and (2) one LIGO detector 
at Livingston, one GEO600 and one VIRGO detector. 
Fig 121 also shows the ratio of the differences in eigenval- 
ues for the two networks. It is clear that the difference in 
eigenvalues for the first network, which consists of nearly 
aligned detectors, is larger than that for the second net- 
work over a large fraction of the sky. 

The statistic Apl(0, </>) solves the two detector para- 
dox by continuously bridging the aligned and misaligned 
detectors cases. From Eq. 021 it is clear that as — Mi 
increases, i.e., as the detectors become more aligned, the 
denominator increases in the negative direction while the 
numerator, sin 203,, stays the same. Hence, 2?/; — > tt 
which implies ip tt/2. From Eq. ^2 the penalty func- 
tion has its lowest value for ip = tt/2. This is the direc- 
tion in the response plane along which all the responses 
in a network get oriented as the detectors are aligned. 
Therefore, with progressive alignment of the detectors 
in a network, the solution ijj is forced towards equal re- 
sponses in all detectors, thus leading smoothly into the 
completely aligned case. In contrast, for the standard 
likelihood statistic Asl(^, (j)) (= ^pl when a = 0), there 
exists a fundamental discontinuity between the solutions 
for aligned and misaligned cases: For the aligned case 
-0 — 7r/2, while ip = (f>x for the misaligned case. 

Finally, the penalty term, F^, depends only the direc- 
tion of the instantaneous network response vector S. It 
does not restrict its length. Hence, the magnitude of S 
at two distinct time instants can be arbitrarily different. 
Thus, even white noise response waveforms are allowed 
as solutions. The only requirement is on the mutual con- 
sistency of detector responses at each instant of time. 



C. Connection with Bayesian analysis 

The penalized likelihood approach, Eq. 021 has a 
straightforward connection to Bayesian ISj analysis. In 
the latter approach, the posterior degree of belief over 
the space of waveforms, or equivalently p and ip (c.f., 
Eq. ref), is obtained using Bayes' law, 

p{p,m)- '^'^''^f^^'^^'\ (50) 

where p(X|p, -0) is simply the likelihood A in Eg. 1111 The 
maximum a posteriori Bayesian estimates for the true p 
and ■(/; would be the ones for which ■i/'IX) is maxi- 
mized. This is equivalent to maximizing the logarithm of 
p{p, ipl'K.) and, hence, the sum of A and logp{p, ip). Com- 
parison with Eg. 1461 shows that 1/F^ acts as a prior. The 
difference is that this prior is not normalized in Eq. 1461 
but is instead multiplied by a. Of course, a can be cho- 
sen to be the integral of 1/F^ over tp thus reducing the 
penalized maximization to a true Bayesian analysis. 



VI. RESULTS 



The PMLR method was implemented in Matlab 
and its performance was quantified with the Monte Carlo 
simulation code previously used in p^ . We give a brief 
recapitulation of the simulation scheme: Binary Black 
Hole merger waveforms from [25^ are injected into short 
stretches of white noise, one waveform per stretch. The 
signals correspond to face-on binaries scattered uniformly 
on the sky with randomly selected polarization angles. 
The signals have a variable matched filtering SNR due to 
their sky position and polarization, but an overall factor 
G is used to scale the amplitudes of the two polarizations 
h+{t) and hx{t). The sky-averaged SNR of the injected 
signals is approximately related to G as SNR ~ 2.36*. 
Further details can be found in [T5l |. 

Consider, first, some case studies to give a qualitative 
feel for the performance of PMLR compared to the SL 
and constraint likelihood methods. Fig ^ compares SL, 
constraint likelihood (hard constraint) and PMLR (with 
a — 1500) for the same signal using the three detector 
network: Hanford, Livingston and GEO. The true sky 
position of the source was 9 = 50°, (j> = 280°. The best 
estimates for sky position (in degrees) were (1) 6 = 76°, 
(j) = 32° for SL and (2) 6* = 48°, = 280° for PMLR and 
constraint likelihood. Note that, compared to SL, the sky 
maps for constraint likelihood and PMLR are less noisy 
and show more contrast between the region in which the 
signal is located and the rest of the sky. 

The simulation code finds out the fraction of all in- 
jected binaries detected by the PMLR (or other) method. 
The detection statistic for all methods is the maximum 
value over the respective skymap (e.g. Fig^. Fig [51 
shows the receiver operating characteristics (ROC) for 
SL, constraint likelihood (both hard and soft constraint) 
and PMLR (with a = 150 and 1500). The measured frac- 
tion of detected signals is plotted against the false alarm 
rate corresponding to the detection threshold used. It is 
clear from the figure that PMLR can vary continuously 
in performance, ranging from SL for low values of a to 
the hard constraint method for high values. 

The simulation code also estimates the error in sky 
position reconstruction. Fig. j^l shows the fraction of the 
injected signals that were detected within 8 degrees of 
their true position, for a range of the scale factor G. The 
PMLR method {a = 1500) has practically the same per- 
formance as the hard constraint method and both per- 
form somewhat better than the soft constraint method. 
The performance of SL, however, is significantly worse 
than all the other methods. This can also be seen quali- 
tatively in the sky maps of the various methods (Fig. ^ 
where the SL sky map always has less "contrast" around 
the true signal location compared to the constraint like- 
lihood and PMLR methods. 
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FIG. 3: The difference in eigenvalues fio and plotted as a function of sky-position angles 9 and </). (All panels display the 
logarithm to base 10 of — Mi-) The top left panel is for a network consisting of three LIGO detectors: two at Hanford and one 
at Livingston. The top right panel is for a network consisting of one LIGO detector at Livingston, one GEO and one VIRGO 
detector. The first network (top left), which is more aligned, shows much larger values of fio — Hi than the second (top right) 
which has more misalignement between the detectors. In the bottom left panel, the part of the sky where — for the first 
network is < 10 times that for the second network has been suppressed. It is clear that this is a smaller fraction of the sky 
compared to the part where — Mi for the first network is > 10 times that of the second network. In the bottom right panel, 
only the part of the sky where /lo — IJ-i for the second network is greater than or equal to that for the first network is shown. 



VII. CONCLUSION 

There exists a serious problem, called the two detec- 
tor paradox, in the formal application of the standard 
likelihood method to the detection of burst GW sig- 
nals. We resolve this problem by imposing a restriction 
on the detector responses allowed as candidate solutions 
for maximizing the likelihood functional. The restriction 
is imposed through a penalizing function that measures 
the variability of the SNR of a candidate solution when 
the corresponding source is displaced slightly on the sky. 
High SNR variability sources lead to very different detec- 
tor responses even when the detectors are closely aligned. 
The standard likelihood method allows such sources as 
perfectly valid candidate solutions. However, such ab- 
normal sources have an inherently reduced observability. 
The proposed method (PMLR) sacrifices the detectabil- 
ity of suuch sources while boosting the sensitivity and 



sky position reconstruction accuracy for normal signals. 

PMLR contains a free parameter, the penalty factor 
a. For a — PMLR simply reduces to the SL method. 
However, method is not very sensitive to the exact value 
of a, and the sky map of the likelihood converges rapidly 
as a is made large. The role of a needs to be investigated 
further. 

Results from numerical simulations show that the per- 
formance of the PMLR method can be tuned, using a, 
to vary over the range covered by the standard likeli- 
hood and the constraint likelihood methods introduced 
in KMRM. For the type of source population used in the 
simulations, a high value of a yields the same perfor- 
mance as the hard constraint method. The ROC curves 
are mildly better compared to standard likelihood but 
the accuracy of source direction reconstruction is signifi- 
cantly better. 

PMLR is also related to the Bayesian approach. Work 
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FIG. 4: Comparison of SL (top left), constraint likelihood (hard constraint) (top right) and PMLR (bottom left) for a = 1500. 
The source is at S = 50°, = 280° and has an amplitude scale factor of G = 3. The network consists of one detector at 
Hanford, one at Livingston and GEO. The detector responses, assuming white Gaussian noise, are shown in the bottom right 
panel. 



in the latter direction has focussed on using priors that 
suppress candidate solutions that are not smooth. (This 
is equivalent to the requirement that the signals have a 
smaller frequency bandwidth than the noise.) However, 
the penalty in PMLR is not based on the smoothness of 
signals in time. It is on signals that lead to widely differ- 
ent responses in nearly aligned detectors. In particular, a 
white noise burst signal is perfectly admissible as a can- 
didate solution in PMLR. If smoothness constraints are 
desired, they can be added on to the PMLR restriction. 

The approach introduced in this paper can be extended 
to include variability of observed SNR as a function of 
source inclination to line of sight and the projected ori- 
entation of the source on the sky. This extension may 
result in a penalizing function with better performance 
than the one used in this paper. 
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FIG. 6: The fraction of signals detected within 8° of their true 
positions. As the signal SNR decreases, all methods lose both 
detection efhciency and incur large errors in direction estima- 
tion. At the same time, as the SNR increases, all methods 
converge to the true signal location in almost all cases. The 
difference between the methods is, therefore, maximum for 
intermediated SNRs. 



